!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate point on free pressure boundary  
! papramenters by MOC.
! but this subroutine usually used follow after a Prantle-Meyer or
! a shock treating
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-18
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : FreeP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   l0       slope of line 34 (streamline)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine FreeP()
   use VariableDef
   implicit none
!                4()
!        3
!             2
!      given variable
!   real*8::x2,y2,Vu2,Vv2,x3,y3,Vu3,Vv3
   integer::icor,i
!      temp variable
   real*8::Vup,Vvp,yp,Vp,thetap,ap,Mp,lp,Qp,Rp,Sp,Tp,Vu30,Vv30,l0
   real*8::x4f,y4f,Vu4f,Vv4f
!      calculate variable
!   real*8::x4,y4,Vu4,Vv4,V4,p4,T4,rou4,M4

   x4f=0.0
   y4f=0.0
   Vu4f=0.0
   Vv4f=0.0
   p4=p0
   
   V4=sqrt(2*g*Rg*Tt*(1-(p4/Pt)**((g-1)/g))/(g-1))

   y2=abs(y2)
   y3=abs(y3)
   Vv2=abs(Vv2)
   Vv3=abs(Vv3)
   Vup=Vu2
   Vvp=Vv2
   Vu30=Vu3
   Vv30=Vv3
   yp=y2

   do i=1,icor
      Vp=sqrt(Vup**2+Vvp**2)
      thetap=atan(Vvp/Vup)
      call Thermo(Vp,p4,rou4,T4,ap,Mp)
      
      lp=tan(thetap+asin(1/Mp))
      Qp=Vup**2-ap**2
      Rp=2*Vup*Vvp-Qp*lp
      Sp=delta*Vvp*ap**2/yp
      
      l0=Vv30/Vu30
      x4=(lp*x2-l0*x3-y2+y3)/(lp-l0)   !location
      y4=y3+l0*(x4-x3)
      
      Tp=Sp*(x4-x2)+Qp*Vu2+Rp*Vv2
      
      Vu4=(Qp*Tp-Rp*sqrt((Qp**2+Rp**2)*V4**2-Tp**2))/(Qp**2+Rp**2)   ! velocity
      Vv4=sqrt(V4**2-Vu4**2)

      if ((abs(x4-x4f) <= el) .or. (abs(y4-y4f) <= el)) then
         if ((abs(Vu4-Vu4f) <= ev) .and. (abs(Vv4-Vv4f) <= ev)) then
            exit
         end if
      end if

!      record as next comparison
      x4f=x4
      y4f=y4
      Vu4f=Vu4
      Vv4f=Vv4
      
!      correct step initial value
      Vup=(Vu2+Vu4)/2
      Vvp=(Vv2+Vv4)/2
      yp=(y2+y4)/2
      Vu30=(Vu3+Vu4)/2
      Vv30=(Vv3+Vv4)/2
   end do
   V4=sqrt(Vu4**2+Vv4**2)
   call Thermo(V4,p4,rou4,T4,ap,M4)   ! calculate paramenters

   if (wallflag==-1) then
      y2=-y2
      y3=-y3
      y4=-y4
      Vv2=-Vv2
      Vv3=-Vv3
      Vv4=-Vv4
   end if

end subroutine FreeP
